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Abstract 

The hydrodynamic forces exerted by a fluid on small isolated rigid spherical particles are 
usually well described by the Maxey-Riley (MR) equation. The most time-consuming 
contribution in the MR equation is the Basset history force which is a well-known problem 
for many-particle simulations in turbulence. In this paper a novel numerical approach is 
proposed for the computation of the Basset history force based on the use of exponential 
functions to approximate the tail of the Basset force kernel. Typically, this approach 
not only decreases the cpu time and memory requirements for the Basset force compu- 
tation by more than an order of magnitude, but also increases the accuracy by an order 
of magnitude. The method has a temporal accuracy of O (At 2 ) which is a substantial 
improvement compared to methods available in the literature. Furthermore, the method 
is partially implicit in order to increase stability of the computation. Traditional meth- 
ods for the calculation of the Basset history force can influence statistical properties of 
the particles in isotropic turbulence, which is due to the error made by approximating 
the Basset force and the limited number of particles that can be tracked with classical 
methods. The new method turns out to provide more reliable statistical data. 
Keywords: Basset history force, numerical approximation, particle laden flow, Maxey- 
Riley equation, isotropic turbulence 



1. Introduction 

The turbulent dispersion of small inertial particles plays an important role in environ- 
mental flows, and in this work we focus on small particles with densities of the same 
order as that of the surrounding fluid. Examples of such particles that may be present 
in well-mixed or in density stratified estuaries arc plankton, algae, aggregates (all with 
densities similar to the fluid density) or resuspended sand from the sea bottom (particle 
densities in this case several times that of the fluid). Particle collisions and the formation 
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of aggregates of marine particles or sediment depend on the details of the small-scale tra- 
jectories of the particles in locally homogeneous and isotropic turbulence. At these scales 
the details of the hydrodynamic forces acting on (light) inertial particles are relevant. 

Maxey and Riley [l[ introduced the equation of motion for small (d p <C rj, with d p the 
particle diameter and rj the Kolmogorov length scale) isolated rigid spherical particles in a 
non-uniform velocity field u(x, t). An important assumption is that the particle Reynolds 
number Re p = c? p |u — u p \/y <C 1, with u p the velocity of the particle and v the kinematic 
viscosity of the fluid. The relative importance of the hydrodynamic forces depends on the 
ratio of particle-to-fluid density and the particle diameter. The computation of all the 
different forces in the Maxey- Riley equation is an expensive time- and memory consuming 
job. Therefore, assumptions are often made regarding the forces that can be neglected in 
the study of particle dispersion. The number of studies underpinning these assumptions, 
however, is rather limited due, for example, to lack of efficient algorithms to take into 
account the effects of the Basset history force with sufficient numerical accuracy. An 
elaborate overview of the work on the different terms in the Maxey-Riley equation and 
their numerical implementation can be found in the paper by Loth Q. 

The term most often neglected is the Basset history force because of its numerical 
complexity. Many recent studies underline the importance of the Basset force compared 
to the other hydrodynamic forces in the Maxey-Riley equation for particle transport in 
turbulent flows, see Refs. [3|, |4|, la, [6| . Moreover, it can affect the motion of a sedimenting 
particle [7( or bed-load sediment transport in open channels, where the Basset force 
becomes extremely important for sand particles [1, Q . It also might alter the particle 



velocity in an oscillating flow field [1 01 ] or modify the trapping of particles in vortices 11 



Fast and accurate com put ation of the Basset force is far from trivial. Although several 



attempts have been made [12|, |13|, |14| , the computation of the Basset force is still far more 
time consuming and less accurate than the computation of the other forces in the MR 
equation. Therefore we present a new method that saves time, memory costs and is more 
accurate. 

The MR equation and the subtlities with regard to the computation of the Basset 
history force are introduced in Section [5J Next, in Section [3] and 21 the new method 
is introduced, where Section [3] focuses on the approximation of the tail of the Basset 
history force and Section 0] on the numerical integration of the Basset history force. 
Thereafter, validation of the method using analytical solutions is discussed in Section [5] 
A simulation of isotropic turbulence, with light inertial particles embedded in the flow, 
has been performed. In Section [5] we compare the results from this simulation with the 
new implementation of the full MR equation with the old version used by van Aartrijk 
and Clercx [B[ . Finally, concluding remarks are given in Section [7] 

2. Particle tracking 

Particle trajectories in a Lagrangian frame of reference satisfy 

with x p the particle position and u p its velocity. According to Maxey and Riley [l[ the 
equation of motion for an isolated rigid spherical particle in a nonuniform velocity field 
u is given by 
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= Fgt + Fp + Fg + Fam + Fb- (2) 



The equation of motion includes time derivatives of the form d/dt taken along the particle 
path and derivatives of the form D/Dt taken along the path of a fluid element. The 
particle mass is given by m p , a is the radius of the particle, p = pv is the dynamic 
viscosity, p and v are the density of the fluid and its kinematic viscosity, m/ is the 
mass of the fluid element with a volume equal to that of the particle and e z is the unit 
vector in the opposite direction of the gravitational force. The forces in the right-hand 
side of this equation denote the Stokes drag, local pressure gradient in the undisturbed 
fluid, gravitational force, added mass force and the Basset history force, respectively. 
The Faxen correction proportional to V 2 u has been included in the Stokes drag, added 
mass and Basset force |lo| . According to Homann et al. [l6| these corrections reproduce 
dominant finite-size effects on velocity and acceleration fluctuations for neutrally buoyant 
particles with diameter up to four times the Kolmogorov scale ry. For the added mass term 
the form described by Auton et al. I?) is used. Moreover, the history force convolution 
function g(t) and its kernel are 

g(i) = ^, f(<) = u-Up + i a 2 V 2 u, = i=. (3) 

Equation ([2]) is valid when a <C rj, but, as mentioned above, the Faxen correction can 
weaken this condition. Furthermore, the particle Reynolds number must be small (Re p <C 
1), as are the velocity gradients around the particle. Finally, the initial velocity of the 
particle and fluid must be equal. The coupled system (fT]) and ^ is in principle suitable 
for integration by any standard method, e.g. the fourth order Runge-Kutta method. 

The Basset history force Fb presents additional challenges. First, the evaluation of the 
Basset force can become extremely time consuming and memory demanding. This is due 
to the fact that every time step an integral must be evaluated over the complete history 
of the particle. Several attempts have been made to solve this problem. Michaelides 
uses a Laplace transform to find a novel way for computing the Basset force. This 
procedure can be used for linear problems, but is not suitable for space dependent velocity 
fields for which the coupled system (fT]) and ([2]) is nonlinear. Another solution is provided 



by Dorgan and Loth [13J and Bombardelli et al. 12]. In these papers the integral is 
evaluated over a finite window from t — t wm until t. This can be represented by a change 
in the kernel of the Basset force. The window kernel is thus defined as 

K fA _j K B (t) for t < t win , (A s 
AwinW ~\ for t>t win . W 

The kernel of the Basset force is decreasing very slowly for t — > oo, thus t W [ n must 
be chosen rather large. For Bombardelli et al. [12] this problem turned out to be 
less important because they used a different kernel, which decreases faster for t — > oo. 
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Although the application of the window kernel saves CPU time, the computation of the 
Basset force is still far more expensive than the evaluation of the other forces in the 
MR equation. It turns out to be approximately 100 to 1000 times more time consuming 
depending on the application. 

A second issue concerns the kernel of the Basset force, which is singular for t 0. 
A standard approach to deal with the singularity of the Basset kernel is to employ 
specific quadrature rules such as the second order Euler-Maclaurin formula (l8l | . Another 
approach is presented by Tatom [l9| who uses a fractional derivative method. This 



approach was tested by Bombardelli efc al. [12| ■ From their results it can be easily shown 



that the integration method with specific quadrature rules has only temporal accuracy 
0(y/At) and that the fractional derivative approach has a temporal accuracy O(At). 
In computations of turbulent flows with particles, other discretization methods involved 
are at least second order. Therefore, it is not sufficient to have a first order integration 
method for the Basset force. 

Our goal is to derive a robust and efficient method for the computation of the Basset 
force that overcomes all the problems mentioned above and to find an approach that is 
suitable for different forms of the kernel. Furthermore, our method must be stable and 
at least second order accurate in time. A third requirement is that it should be less time 
consuming and memory demanding than previous methods. 



3. Approximation of the tail of the Basset force 

To get a better understanding of the Basset force we will first show that the contribution 
of this force is finite at any given time. To do this, some restrictions on f(t) and g(t) = 
^f(t) should be made. First, f(t) must be a continuous function and its derivative must 
exist almost everywhere. Further, f(t) and g(t) must be in the L°° space with norm B\ 
and B2, respectively. The restrictions on f(t) and g(t) are thus: 

feC°, \\f\\ 00 = B 1 , ||g|U = S 2) (5) 

where || • |joo is defined as: 

Hfllco = inf{C > : |f(t)| < C almost everywhere}, (6) 

and I ■ I is the usual length of the vector. We assume that for particles in (turbulent) 
flows with f(t) = u — u p + ia 2 V 2 u these conditions are satisfied as both the flow field 
and its Laplacian satisfy these conditions. With the conditions in ([5]) it is possible to 
find an upper bound for F^. The integral is split into two parts, in order to control both 
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the singularity in the Basset kernel and the tail of the integral. This yields 
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Here cb = d>a 2 py/ : Kv is introduced for convenience. We now consider the window kernel 
for calculation of the Basset force Fs-win ■ In the limit of i w in — > oo the difference between 
Fb and Fs-win must vanish. Using integration by parts, one can derive 
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The error made by using the window kernel instead of the Basset kernel is indeed becom- 
ing negligibly small for i w ; n — > oo. Unfortunately, this convergence is very slow, implying 

that twin 

must be very large, and a better approach for the computation of the Basset 
force must be found. This is done by introducing a new kernel with a modified tail, in 
short the modified Basset kernel K mo d(t) 7 as follows 



K 



mod 



(*) 



K B (t) for t<t win 

K tai l(t) for t > t w in 



lim K tail (t) = 0. 

t—>00 



(9) 



This new kernel also implies a modified history force denoted by FB-mod- For now K ta n(t) 
is not yet defined but must be chosen such as to approximate the Basset kernel as close 
as possible. Using integration by parts in the last step, the upper bound for the error 
induced by the modified Basset force FB- mo d becomes: 
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-Kb (twin) — -Ktail(twin) 
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d(K B - K tail )(t) 



dt 



dtj (.10) 



As the upper bound in relation fj 10[) depends on i W i„, it turns out to be beneficial to 
rescale the time and kernel as follows: 



-Ktail(t) 



K~q win ) ^win 

Applying the same scaling to Ks(t) = we find 



(11) 
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-f^B^win) 

Note that this cannot be done for a general kernel. Eq. (fit))) can now be reformulated as 
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(13) 



When comparing ([5J and (fT3|) one can see that a good approximation K ta ,ii(t) of the tail 
reduces the error in (|13p significantly in comparison with ([5]). 

In order to find a good approximation iftail(^) we start with (|10[) . The right hand side 
of (jlOp can be minimized and thereby minimizing the error in FB_ m od ■ When determining 
-Ktaii(0 it is important that computation time is kept low. In order to achieve this, 
exponential functions are used because they can be implemented in a recursive way as 
explained later on. At first we start with one exponential function as follows, 



K ta ,n(t) = aexp(-ta) . 



(14) 



Here a and b are two positive constants. As a first guess we require that -Ktaii(iwin) 
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^B(^win) in order to determine a and b. In this way 



K mo d(t), defined in ([9]), is continuously differentiable. Doing this results in 



#tail(f) = \/r^-exp ( - — 

twin \ ^^wi 



(15) 



Fig. Q] shows several kernels, where the modified Basset kernel is given by (|15l) . The 
error by applying the modified Basset kernel is obviously smaller compared to the error 
for the window method. In order to minimize the error even more, multiple exponential 
functions can be used. Relation (fl~5l) provides an ansatz for the choice of a and b. Thus 
we write K t nn(t) as 



K tai i(t) 



E 



exp 



t 

2t, 



(16) 



with a,i and t{ positive constants. The functions Kiit) satisfy the following properties: 
Ki{ti) = K B (U) and gj-Ki(ii) = ^K B {U). Combining (JTTJ) and ([TB]). we obtain the 
following dimensionless representation for the tail: 



K tai i(t) 
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Figure 1: Basset kernel (solid line), window kernel (dots) and the modified Basset kernel (dashed line) 
for twin = 2. 

The coefficients a.; and U should be chosen in such a way that the upper bound in (|13[) 
is minimized. However, Newton iteration will not work for this problem, and instead we 
consider the expression 



which provides a good indication for the optimal values of ai and tj. In (|18[) an extra 
multiplication with t is introduced to correct for the change in norm. After minimizing 
the expression in (|18[) . we can verify whether the error in fj 13[) is of the same order. Since 
Ki(ti) = Kb(U) and ^-Kifa) = 4?Kb (£»)) the function Ki(t) approximate K&(t) very 
well around U. The kernel Kb must be approximated over a large range of t- values and 
as a consequence i, must also have a large range. Furthermore, Kb is changing slowly 
for large t so the small t, must be close to each other whereas the large t, can be far 
apart. The approach for finding <Zj and ti is thus the following. First, make a reasonable 
choice for ti, and second, calculate at by minimizing (I18[) . Finally, determine the term 
between brackets from (| 1 3j) . Another slightly different set of ^-values can be chosen to 
see if a better approximation can be made. In Tabled] the result is shown for m = 10. 
Here one can see that some values of U are smaller than 1. This is surprising because the 
kernel Kb is not being approximated below t = 1. When tuning the ij-values we found, 
however, that this improves the approximation. 

From Fig. [5] it can be seen that -Ktaii approximates Kb relatively well over a wide 
range of t. From Fig. [3] one can see that the error decays for large t (note the huge range 
of t in both figures) . 

Using (jTTJ) in combination with Table Q] for K ta i\(t) the part between brackets in (fT5|) 




(18) 
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Table 1: Coefficients a; and ti in -KtailM with m = 10 



u 


a, 


0.1 


0.23477481312586 


0.3 


0.28549576238194 


1 


0.28479416718255 


3 


0.26149775537574 


10 


0.32056200511938 


40 


0.35354490689146 


190 


0.39635904496921 


1000 


0.42253908596514 


6500 


0.48317384225265 


50000 


0.63661146557001 
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can be calculated 

1 - ^tail(l) 

Comparing this result with the window method ([8j a factor of more than 200 is gained 
in accuracy. When keeping the same accuracy but changing the window, i w j n can be 
decreased by a factor of 200 2 = 40000. 



d(K B - K tan )(t) 



dt 



dt rj 9.5 • 10 



-3 



(19) 



4. Numerical approximation 

In this section the numerical integration is discussed. First, the integration of the window 
and tail kernels are elaborated. Second, the overall numerical scheme for solving Eq. (JXJ) 
and ([2]) is explained. 

The integration of the Basset force with the modified kernel ([9]) and (fT6|) is split into 
two parts, the window kernel and the tail of the kernel as follows, 



B-mod 



(t) = c B / K mod (t - r)g(r)dr 

J — oo 

= c B / K B (t - T)g(r)dT + c B / K tai \(t - r)g(r)dr 

^t-t„i„ J-oo 
= F B -wm(i)+F B -tail(i)- (20) 

In the following, methods are described for the calculation of F B _ win and FB-taii- 

First, we consider the Basset force due to the window kernel F B _ win . The kernel of 
the Basset force is singular for t — > which impedes use of the ordinary trapezoidal 
rule. In order to deal with the singularity we introduce an alternative, trapezoidal-based 
method, referred to as the TB-method. The idea is as follows. The trapezoidal rule 
is based on linear interpolation of the integrand on each subinterval. In our approach 
g(t) is approximated by its linear interpolant Pi(i), and subsequently the integration of 
K B (t — t)Pi(t) is done exactly. For the numerical implementation we start with the 
discretization of the interval [t — twin, t], given by r„ = t — nAt, n = 0, 1, 2, • • • , N with 
At = t w i n /N. Now the integral can be split as 

F B -winW = c B V f"" ~^=dr . (21) 

The next step is to approximate g(r) by its linear interpolant on each subinterval, which 
yields 

Ffi-winW ~ c B > / n dr, (22) 

where g n = g(r„). After the change of variable r' = t — r this integral can be evaluated 
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and the following resul10 is obtained: 
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(23) 



From the result above one can see that three inner products must be calculated each time 
step, one inner product for each spatial dimension. One vector contains all the values 
g n which must be shifted by one index each time step. The other vector containing the 
coefficients in (|23p is calculated once at the start of the computation. In this way the 
computational time is kept minimal. The part with g will be treated in a different way 
as explained later on in order to improve stability. 

Next, the numerical integration of the tail of the Basset force is discussed. The idea is 
to find a recursive formulation in order to minimize computation efforts. Using expression 
(HSJ) for if tail, F B -taU becomes: 



' B-tail 



(*) 
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CljCB 



(24) 



Here, Fj represents the contribution of the i-th exponential function. Now Fj is split 
into two parts, as follows. 



Fi(t) = c B 



Ki(t - r)g(r)dr + c B 



-At 



t-twin-At 



K t (t - T)g{ T )dT 

= F i . di (t)+F;_ re (i) , (25) 

where we have to compute F;_di directly and where Fj_ ro can be computed recursively. 
For Fi_di the same procedure is followed as with the window kernel. Using this procedure 
the following result^ can be obtained: 
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Sn 
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. (26) 



where (f(z) = (e z — l)/z = 1 + \z + ^z 2 + O (z 3 ). Finally, Fi_ re can be easily calculated 
using the value of Fj at the previous time step: 

r-t-twin-At 



Fi_ re (t) 



CB 

exp 
exp 



exp 
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2f, 
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t-At-r 



■ exp 



2t, 



g(r)dr 



Fi(t- At) . 



(27) 



1 This formulation is preferred to avoid loss of significant digits in the computation of Fb-- 
2 Note that in equation (1 26 I t Taylor series must be used for ip ^— ^jM when At <C tj. 
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In this last part the overall numerical scheme is discussed. To solve equation (JTJ) and 
([2|) numerically the second-order Adams-Bashforth (AB2) method is implemented. For 
a differential equation 4? = h(t, y) the scheme reads y„ +1 = y n + ^ (3h" — h" -1 ), 
where h™ = h(i™,y"). Equation ([T]) can be directly integrated with this scheme but for 
equation @ some modifications are needed. In order to have a stable scheme, the 
term in the added mass force is treated in an implicit way instead of explicit. Moreover, 
it turned out that the AB2-method has poor stability properties for the calculation of 
the Basset force using the window method. Extremely small time steps must be taken in 
order to have a stable solution. An alternative method circumventing stability problems 
is to bring a part of the Basset force (the contribution evaluated at t) to the left 
hand side. Eq. © is then reformulated as 

v + \ m f + ^bVA^) ^ = F St + F P + F G + F' AM + F B , (28) 

with F' AM = \m f + ^a 2 ^(V 2 u)) and F' B = F B - \c^\fKt^-. In this way the 
Basset force becomes partially implicit instead of completely explicit. Finally, as only 
the time derivative along the particle path ^ is available, the time derivative along the 
path of a fluid element ^ is computed according to 

Du cHi du du du du du du 

m = m +u ^ = ^ + Up ' 3 d^ + {Uj ~ Up ' j) d^ = dF + {U] ~ Up ' 3) d^ ' (29) 



5. Validation of the Basset force integration 

In this section four test cases are presented in order to validate the methods for the inte- 
gration of the Basset force. The first example tests the trapezoidal-based (TB) method 
and compares the results with the semi-derivative (SD) approach by Bombardelli et aJ. 
[l^ . Example 2 and 3 test the the overall numerical scheme. Here both stability and 
convergence are tested for the explicit and the partially implicit TB-mcthod. Finally, 
example 4 shows the efficiency of the Basset force using the tail kernel. 

Example 1: Basset integral for a given convolution function. 

In order to demonstrate the advantages of the TB-method, the convergence of this 
method is compared with the SD-approach of Bombardelli et si. To that end 

the arbitrary test function g(r) = cost is used. The exact Basset integral is given by 

ft cosr 

F B (t) = c B I . dr = 2c B / cos(t- er 2 )dcr 
Jo Vt-T Jo 

= c B \^ (^C\ v / 2l/^)cost + S(y / 2t/^)smtj , (30) 

with a = \A — r and C(t) and S(t) the Fresnel cosine and sine functions [2(|, respectively. 

The Basset integral F B was evaluated at t — 50ir with different numbers of points 
N uniformly distributed in the interval [0, 507r]. The results for both the SD-approach 
and the TB-method are presented in Table [2] Here, it can be seen that the error of 
the TB-method is substantially smaller than that of the SD-approach. When increasing 
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the number of points N it can be seen that the TB-mcthod is second-order accurate 
in time (in agreement with analysis that can be done by using Taylor series), whereas 
the SD-approach is first-order accurate in time. More methods have been compared by 



Bombardelli et al. [12j but these methods have even lower order of convergence than the 
SD-approach. 



Table 2: Relative error and order of convergence for the Basset integral, for the SD-approach 1 1 2l and 
the TB-method. 





Relative error 


Order 


Relative error 


Order 


points N 


SD 


SD 


TB 


TB 


81 


4.03- 10" 1 




1.34- lO" 1 




243 


1.37- 1CT 1 


1.0 


2.54-10- 2 


1.5 


729 


4.66- 1(T 2 


1.0 


3.29-10" 3 


1.9 


2,187 


1.56- 1(T 2 


1.0 


3.93- 10~ 4 


1.9 


6,561 


5.22 • 1(T 3 


1.0 


4.54- 10~ 5 


2.0 


19,683 


1.74- 1(T 3 


1.0 


5.15 • 10~ 6 


2.0 


59,049 


5.80 • 10~ 4 


1.0 


5.80- 10- 7 


2.0 


177,147 


1.93- 10~ 4 


1.0 


6.49 • 10~ 8 


2.0 


531,441 


6.45 • 10~ 5 


1.0 


7.24- 10~ 9 


2.0 


1,594,323 


2.15 • 10~ 5 


1.0 


8.06- 10~ 10 


2.0 



Example 2: Space- dependent steady velocity field. 

In order to test the overall numerical scheme for the computation of particle trajectories 
we have implemented a particular space-dependent steady velocity field. The particle 
trajectory is a circle and given by (x(t),y(t)) — (r cos cot, — r smut), where r and u denote 
the radius and the angular velocity, respectively. The velocity field and its derivation is 
given in appendix A. For the test case, exactly one revolution is simulated, from t = 
until t = 2ir. In order to test the stability of the overall scheme two different approaches 
have been tested. One with the completely explicit time integration procedure for the 
Basset force and the other with the partially implicit procedure, see Section 2] For both 
the implicit and explicit method the Basset force is computed with the TB-method and 
show second-order convergence in At. The relative error is computed with x p (27r). The 
results are presented in Table [3] and clearly indicate that the explicit scheme is very 
unstable when the number of time steps is smaller than 256. The partially implicit 
scheme remains stable even with the number of time steps as small as 16. 

Example 3: Time- dependent velocity field. 

The trajectory of a spherical particle in an arbitrary time-dependent velocity field can 
rather straightforwardly be computed as long as the velocity field is smooth enough. 
The derivation of the particle trajectory uses Laplace transforms and the analytical 
procedure is given in appendix B. The overall numerical scheme is tested by computing 
the trajectory of a particle in the following one-dimensional, time-dependent velocity 
field 

, . (m v — mAq „ . , 

u(t) = s — cos 2t . (31 
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Table 3: Relative error and order of convergence for the overall numerical scheme, tested for the trajectory 
of a small particle in a space dependent steady velocity field. 



number of 


Relative error 


Order 


Relative error 


Order 


time steps 


explicit 


explicit 


implicit 


implicit 


16 


unstable 




3.63- IO" 1 




32 


unstable 




8.32- 10~ 2 


2.1 


64 


unstable 




2.09- 10~ 2 


2.0 


128 


unstable 




5.28- 10~ 3 


2.0 


256 


4.80- ia 2 




1.33- IO" 3 


2.0 


512 


3.05- 10~ 4 


7.3 


3.33- 10~ 4 


2.0 


1024 


7.68 • 10~ 5 


2.0 


8.34- 10~ 5 


2.0 


2048 


1.93- 10~ 5 


2.0 


2.09- 10~ 5 


2.0 


4096 


4.84- 10~ 6 


2.0 


5.21 • 10~ 6 


2.0 



The total force on the particle is zero at t = 0, i.e, Fgt and Fg are in balance. In 
order to compute the Basset force the implicit TB-method is used. The integration is 
carried out from t = until t = 2tt. The relative error is computed for u p (2tt) and is 
presented in Table 21 where once again second-order time accuracy is confirmed. From 
these test cases, using both a time-dependent and a space-dependent velocity field for the 
computation of particle trajectories, we can conclude that the (partially implicit) TB- 
method is stable and second-order accurate in time, and conjecture that this remains the 
case for particles in arbitrary time- and space- varying flow fields. 

Table 4: Relative error and order of convergence for the overall numerical scheme, for the velocity field 

(31) ■ 



time steps Relative error Order 



16 


9.96 


10" 2 




32 


2.38 


io- 2 


2.1 


64 


5.57 


10~ 3 


2.1 


128 


1.31 


10~ 3 


2.1 


256 


3.13 


io- 4 


2.1 


512 


7.56 


10~ 5 


2.0 


1024 


1.84 


10~ 5 


2.0 


2048 


4.53 


IO" 6 


2.0 


4096 


1.12 


10" 6 


2.0 



Example 4-' Computational efficiency due to modified kernel integration. 
In this example the computational savings when using the modified tail kernel, given in 
© and (fTo) . is investigated based on analysis of the number of flops per time step, per 
particle and per space dimension. For the window kernel this is N + 1 flops because only 
one vector dot product is calculated. For each exponential function three extra flops are 
needed. To see how efficient the tail kernel works the upper bound (fl"3"|) for the error is 
plotted as a function of the computation time, Fig01 Different numbers (indicated by m) 
of exponential functions are taken into account. The results are plotted in Fig 01 Here it 
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can be seen that the best choice for m depends on the particular situation. Furthermore, 
the results show a significant saving in computation time. This can easily be a factor 
of 100 or more. When looking to the memory requirements the results are even better. 
For the window method as many memory locations as flops are needed whereas each 
exponential function only takes one memory location instead of 3 flops. So using the tail 
kernel not only saves time but also memory. 

Overall, the use of the tail kernel reduces the computational costs of the Basset force 
by more than an order of magnitude, whereas the memory requirement is even reduced 
more. Furthermore, the error is reduced by more than an order of magnitude. The 
question remains, of course, whether the computational savings directly result in faster 
simulations. This depends on the remaining part of the simulation. Although the other 
force contributions in f|28[) can be calculated much faster than the Basset force this does 
not have to hold for the interpolation of the velocities in a turbulence simulation. The 
velocity of the flow field is only computed at the grid points and an interpolation must 
be carried out to compute the velocity at the particle position. This may be very time 
consuming and it can become the new bottleneck. The reduction of CPU-time might 
then not be as big as expected but it remains significant. Additionally, the decrease in 
memory requirement may become essential when increasing the number of particles in 
turbulence simulations. 

10° 



10- 2 
I icr 3 



10- 4 
io- 5 



10° 10 1 10 2 10 3 10 4 10 5 

number of flops 

Figure 4: Upper bound for the error in the approximation of the Basset force as a function of the number 
of flops for different number of exponential functions (indicated by m). 




6. Light particles in isotropic turbulence 

In this section a brief statistical analysis of velocities of particles, released in an isotropic 
turbulent flow, is provided. The isotropic turbulence simulation is performed by means 
of direct numerical simulations. The numerical code consist of two parts. First the 
Navier-Stokes equations with the Boussines q appro ximation are solved on a triple pe- 
riodic domain using a pseudo-spectral code [2lj, [22[ (Eulerian approach). Second, the 
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particle trajectories are obtained by the Lagrangian approach as explained in the previ- 
ous sections. The simulation is performed on a 128 3 grid. The number of (light) particles 
is 20,000 and the particle-to-fluid density ratio Pp/pf — 4 (thus all hydrodynamic forces 
in the MR equation are relevant, see Refs. 0, The integral-scale Reynolds number 
is Re = UL/v = 1333, with U the typical root-mean-square velocity and L the integral 
length scale. The Stokes number St is typically in the range 0.1 < St < 1.0 @ and 
particles are tracked for a period of approximately two eddy turnover times. 

Two simulations have been carried out under exactly the same flow conditions and 
particle tracking is either based on the classical approach (window method) or on the 
novel integration method (exponential method) for the Basset kernel. In the first simu- 
lation only the window kernel Q has been used, where the number of time steps in the 
window is n — 500. The other one uses the modified window kernel, given in ([9]) and 
(I16|) . In this case only five time steps are taken into account in the window, so n = 5. 
For the tail of the Basset kernel the number of exponential functions m = 10. 




Figure 5: Autocorrelation of the particle velocity 
Up. The solid line represents the result from the 
exponential method and the dots those from the 
window method. 



Figure 6: Energy spectrum of the particle velocity. 
The graph of the window method (dashed line) is 
shifted downward with respect to the spectrum 
from the exponential method (solid line) by a fac- 
tor of 100 for clarity. 



In order to study a particle trajectory we start with considering the energy spectrum of 
the particle. To obtain the energy spectrum, we first need to calculate the autocorrelation 
R(t) of the velocity, which is defined by 

_ (u p (t)u p (t + t)) 

r(t) - Mm ■ (32) 

Here, (•) denotes the average in over the different particles. The particles are embedded 
in a homogeneous isotropic turbulent flow and no gravitation is applied. Therefore, we 
are allowed to average over the components of the velocity vector of all particles. No 
time averaging has been applied for the present velocity data as this run covers only one 
or two eddy turnover times. The results for the autocorrelation of the velocity are shown 
in FigJS] and we see that the results for both the window method and the exponential 
method are comparable. The energy spectrum obtained from the particle velocities can 
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be calculated by taking the cosine transform of the autocorrelation function, and is shown 
in FigJHl Although the results are similar we are interested in possible differences between 
the two spectra. If these differences have an overall trend this would mean that statistical 
properties can be influenced by the different methods of evaluating the Basset kernel. 
However, to observe any error in the evaluation of the Basset force kernel the differences 
should be larger than the statistical noise. 

A starting point for an analytical evaluation of possible differences between the win- 
dow method and the exponential method consists of the response of a single particle 
in a uniform oscillating flow field. We are therefore interested in the periodic solu- 
tion Up of a spherical particle responding to an oscillating velocity field u = coswi (or 
u = TZ[exp(iu>t)], with i the imaginary unit and 1Z denoting the real part of this ex- 
pression). The particle velocity can then be expressed as u p = 1Z[V exp(wji)] with V 
a complex amplitude, which is dependent on the method chosen to evaluate the Basset 
force kernel. For the window method and the exponential method we introduce V w i n and 
V cxp , respectively. For V exp the exact solution V ex is used since the error of the exponen- 
tial method is assumed to be negligibly small, see also Fig. [3] In general, |V w i n | 7^ |T4 X | 
which means that some frequencies are suppressed with the window method while oth- 
ers may be amplified. This should become visible in the energy spectrum of particle 
velocities. 

In order to find V^n Eq. @ should be solved for u — lZ[exp(iu>t)] and u p = 
7?.[V r exp(ia;i)], resulting in the following intcgro-differential equation: 

iLum p V win = 67ra/i (1 - V win ) + y ™/ (3 - V win ) 

r-t e -iuj(t-T) 

+ iuc B (l - V win ) dr . (33) 

Jt-t win y/t-T 

Here, we used the fact that the velocity field is uniform, one dimensional and that no 
gravity is applied. Applying the change in variables a = \J(t — t)lj, allows us to find an 
expression for V w - m i.e., 

V win = 1 + n (34) 

6ira^ + (5m/ + m p ) iui + cbv 2umQ(y'2t win um) 

where Q(t) = S(t) + iC(t), with C(t) and S(t) the Fresnel cosine and sine functions, 
respectively [20I ] . V ex can now be found by taking V cx = lim twin _ i . 00 Vwin which results in 

y ox = i + ( m / ~ m p) luJ . (35) 

6Tra/i + {j^nif + rri p ) iui + cb-v/^(1 + i) 

Inspection of the energy spectrum displayed in Fig.[5]reveals that the noise becomes more 
important for increasing uj. The effects from the different methods for the computation 
of the Basset force kernel turned out to be most important for u) > 1. Unfortunately, 
the noise in the energy spectrum is already larger than predicted for the differences 
between the window and exponential method. One way of decreasing the error would be 
averaging over time, but with the limited number of eddy turnover times in the present 
simulation this is not feasible. However, an alternative approach exists in comparing the 
autocorrelation of the particle acceleration. 
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Figure 7: Autocorrelation of the particle acceler- 
ation a p = d\ip/dt. The solid line represents the 
result from the exponential method and the dots 
those from the window method. 



Figure 8: Spectrum of the particle acceleration. 
The graph of the window method (dashed line) 
is shifted downward with respect to the spectrum 
from the exponential method (solid line) by a fac- 
tor of 10 for clarity. 




Figure 9: The theoretical ratio ( \V\ J (dashed line) compared with a similar ratio of the particle 
acceleration spectra (solid line). 



The autocorrelation of the particle acceleration is plotted in Fig. [JJ Here, the typical 
time scale is much shorter than that of the particle velocity, therefore it is possible 
to also average over time. The spectrum is calculated by taking the cosine transform 
of the autocorrelation acceleration function and is displayed in Fig. [5] Because the 
particle acceleration is used instead of the particle velocity, higher frequenties (shorter 
time scales) become more important. In this way deteriorating influence of the noise 
on the spectrum is shifted to higher frequenties. Nevertheless, the computed spectrum 
should still be affected by the method that is chosen for the evaluation of the Basset 
force kernel. In order to observe the differences the best approach is to plot the ratio 



17 



of the two spectra as function of frequency. When no essential difference exists between 
the window and exponential method their ratio would be equal to one with some noise 
added to it. However, the window method suppresses some frequency components while 
others are amplified, so the deviation from one is a measure for the error in the window 
method. In Fig. [H] the ratio of both spectra is shown in combination with the theoretical 



the ratio obtained from the simulation, including the local maxima and minima quite 
well, provided the frequency is not too high. For higher frequencies the noise becomes 
larger but the theoretical and computational ratios still seem to have the same trend. 
The novel exponential method to evaluate the Basset force kernel might be considered 
as an excellent and efficient method for tracking of many particles in turbulent flows. 

7. Conclusions 

We have introduced a novel method for the evaluation of the Basset force kernel and 
analysed several aspects of its implementation. The tail of the Basset force kernel is 
approximated by exponential functions. The contribution of these exponential functions 
can be calculated in a recursive way which makes it very efficient. Typically the use of 
the tail kernel reduces the computational costs of the Basset force by more than an order 
of magnitude, whereas the memory requirement is reduced even more. Furthermore, the 
error in the tail of the Basset force is also reduced by more than an order of magnitude 
in comparison with the traditional window method. 

A trapezoidal-based method is developed in order to deal with the singularity of the 
Basset force. This method has a temporal accuracy of O (At 2 ) where other methods only 
have a temporal accuracy of O (At) or lower. This method is made partially implicit in 
order to make it more stable. 

The method has been implemented in a tracking algorithm for (light) inertial parti- 
cles in turbulent flows. The isotropic turbulence simulation shows that the error made 
by the window method can influence statistics on the particle trajectories. This has been 
illustrated with the velocity and acceleration spectra. Therefore, the novel exponential 
method is preferred over the classical window method. Because the new implementa- 
tion is much faster than the classical one, more particles can be taken into account in 
simulations, which opens possibilities for further research. 
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Appendix A. Flow field for circular particle trajectories 

In this appendix the space dependent velocity field is derived that allows a circular 
particle trajectory as solution of the MR equation. Suppose the particle trajectory and 
velocity is given by 




From Fig. |H] it can be seen that the theoretical ratio predicts 



x p = rTZ [' 



e e 




te e 



(Appendix A.l) 



with e = e x — ie y . For the flow velocity field we are looking for solutions of ([2]) of the 
form 

u = —TZ [s(x + iy)e] + 2aze z , (Appendix A. 2) 

with s = a + i/3 a complex constant. For (3 = the velocity field represents a sink flow 
u = (—ax, —ay, 2az) and for a = it represents solid body rotation: u = (f3y, — fix, 0). 
A spherical particle released in the plane z = will remain there due to the symmetry 
of the flow. Substituting ( |Appcndix A.l[ ) and ( |Appcndix A.2[ ) (assuming z — 0) in 
equation and taking into account that no gravity is applied, the Faxen corrections 
are for a linear velocity field and = 7Z [s 2 (x + iy)e \ , yields the following quadratic 
relation for s: 

— uj 2 (^rn p + 2™/^ = 67ra/^(— s + iuj) + —rrif-s 2 + CB \j ( w * s ) ((Appendix A. 3) 

There are two solutions for s, but one of the solutions of s results in an unphysical particle 
trajectory and is therefore discarded. 



Appendix B. Time dependent velocity field 

In this appendix the particle trajectory is derived given the uniform, time dependent 
velocity field (|3"Tj) . The particle is released with an initial velocity u p (0). For a uniform 
velocity field Eq. ([2]) can be simplified to 

1 A dw , du 

i p + -mf 1 — = 67ra/iw + (m/ - m P )-^ - ( m P - mf)ge z 

dw(r) 

+ c B / Kb (t — t) — r-^dr, (Appendix B.l) 
Jo dr 

where w = u — u p . The velocity field u will be expanded in a Fourier series, u(t) = 
Sj^L-oo u ne mwt . The Laplace transform of w is given by W(s) = / °° e~ st -w(t)dt, and 
the Laplace transform of equation ( | Appendix B.l[ ) reads 

m p + im/ + c b^) (sW - w(0)) = 67ra/iW - ^—^lge z 

oo 

, . \ - muj 
I taf-mJ > u„ ; . (Appendix B.2) 



n— — oo 



Using spitting in partial fractions this yields for W(s) 

W(s) = — ( - 7 = n ^ ; ; - ) ( w(0) g-^ J . 9 e z ) + g— J ^e 2 

OO 

+ E 



s Ws + c_ A /s + c i 



n= — oo 



D7ra/i / D7ra/i s 
1 \/ inuj(c + + c_) 



(c + + V inuS) [c- + v inw)(s — inuS) (c^_ — inui)(c?_ — incu)y/s(y/s + v iw) 



*s \(c].-!Jiu)(\/s + C|) (c~ — inuj)(^/s + c_) 
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, (Appendix B.3) 



with c+ , c_ , c and c„ constants given by 



c± = 
2m p + mf 



cb^/tt ± ^/c|7r — 12na^(2m p + mj) 



2-\/ c B 7r — 127ra/i(2m p + m/) 
Transformation back to physical space results in 



2m p + mf 
C = »^(^- ^/)> ppendix b.4) 



oo 

+ £ < 



w(0) | ^e, 



67ra/i 



1 



'»(c + + c_) 



(c+ + Vlnw) (c_ + v/inw) ( c + _ inu){c 2 _ - into) ' 

%l> (c + Vtj - -j——. — tp (c-Vi^ | , (Appendix B.5) 



+ - 



with vp(z) = exp (z 2 ) erfc (2;). 
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